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Abstract 

The  purpose  of  this  thesis  is  to  examine  the  secular  effects  of  the  J2  oblateness 
perturbation  on  close  proximity  satellites.  The  main  objective  is  to  analyze  the  deputy’s 
position  and  velocity  with  respect  to  the  chief  and  adjust  the  initial  conditions  of  the 
deputy  in  an  attempt  to  minimize  the  secular  effects  of  J2  perturbations.  Previous  work 
has  provided  a  method  of  obtaining  a  closed  fonn  solution  for  J2  invariance  with  co- 
planar  orbits.  Therefore,  this  work  will  primarily  consider  deputy  orbits  that  experience 
motion  outside  of  the  chiefs  orbital  plane. 

Upon  detennining  the  required  initial  conditions,  the  invariance  will  be  verified 
through  numerical  integration.  The  method  will  be  considered  successful  when  it  is  able 
to  reduce  secular  effects  to  near  numerical  tolerances. 
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MINIMIZING  SECULAR  J2  PERTURBATION  EFFECTS  ON  SATELLITE 


FORMATIONS 


I.  Introduction 


Objective 

The  concept  of  satellite  fonnations  opens  the  possibility  for  multiple  applications. 
One  important  example  is  distributed  aperture  surveillance.  Multiple  sensors  at  different 
locations  may  observe  an  object,  and  the  multiple  vantage  points  allow  the  images  to  be 
refined.  The  result  is  a  finer  resolution  than  the  sensors  could  provide  individually. 
Another  advantage  of  a  cluster  of  small  satellites  would  be  the  distribution  of  both 
capabilities  and  redundant  systems.  This  would  increase  survivability  and  reliability  by 
distributing  redundant,  mission  critical  systems  amongst  several  smaller,  individual 
satellites. 

There  could  also  be  more  flexibility  when  forced  to  reassign  assets.  If  the 
formation  has  mission  capable  redundancies  distributed  amongst  several  smaller 
satellites,  operators  would  have  the  option  of  reassigning  redundant  systems  while 
keeping  the  primary  assets  on  the  current  task.  If  the  formation  had  a  distribution  of 
capabilities  amongst  several  satellites,  this  could  also  allow  flexibility  in  the  case  of  rapid 
reassignment.  In  this  situation  would  be  possible  to  deploy  mission  critical  assets 
quickly.  This  would  allow  these  assets  to  accomplish  more  time  sensitive  tasks. 
Meanwhile,  if  the  other  assets  in  the  cluster  are  less  time  sensitive  a  smaller,  more  fuel 
efficient  maneuver  can  be  made.  This  would  allow  for  the  completion  of  mission 
objectives  without  consuming  more  fuel  than  necessary. 
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While  there  are  many  benefits  of  satellite  formations,  there  are  some 
disadvantages,  as  well.  One  disadvantage  is  the  amount  of  fuel  required  to  maintain  the 
formation.  In  the  presence  of  perturbation  accelerations,  formations  require  significant 
fuel  consumption  in  order  to  maintain  the  desired  relative  trajectories.  For  a  Low  Earth 
Orbit  (LEO),  the  primary  perturbations  are  due  to  the  atmospheric  drag,  solar  pressure, 
and  the  earth’s  oblateness.  While  drag  forces  are  an  important  factor  for  individual 
satellite  dynamics,  the  satellites  would  be  affected  similarly  if  the  satellites  were  identical 
and  maintained  the  same  orientation.  Therefore,  the  relative  position  and  velocity  will 
not  be  greatly  affected.  The  perturbation  acceleration  due  to  the  oblateness  of  the  earth, 
commonly  referred  to  as  the  J2  perturbation,  has  the  greatest  effect  on  the  relative  motion 
and,  effectively,  the  rate  of  fuel  consumption. 

It  has  been  proposed  that  by  adjusting  the  initial  conditions  of  the  formation,  the 
perturbing  accelerations  may  be  reduced.  Because  the  primary  perturbation  acceleration 
for  LEO  orbits  is  the  L  oblateness  perturbation,  the  objective  of  this  thesis  is  to  determine 
a  method  of  adjusting  the  initial  conditions  to  minimize  the  secular  effects  of  the  L 
perturbation. 

Approach 

This  thesis  will  attempt  to  minimize  the  secular  effects  of  the  J2  perturbation 
accelerations  while  keeping  in  mind  the  desired  relative  geometry.  The  problem  will  be 
presented  in  the  fonn  of  a  fixed  chief  orbit  and  desired  Relative  Orbital  Elements,  and  the 
solution  will  include  a  process  of  making  changes  to  the  desired  set  of  Relative  Orbital 
Elements  in  order  to  achieve  near  J2  invariance. 

In  order  to  accomplish  the  objectives  stated  above,  the  approach  taken  in  this 

research  was  to  first  describe  the  secular  effects  of  the  h  perturbation  on  the  formation 
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parameters  as  a  function  of  the  relative  geometry  and  the  chiefs  Classical  Orbital 
Elements  (COE).  The  parameters  used  to  describe  the  relative  geometry  were  Relative 
Orbital  Elements  (ROE),  an  intuitive  set  of  parameters  that  describe  the  size,  shape,  and 
orientation  of  the  motion  relative  to  the  chief,  in  a  manner  analogous  to  how  the  COE 
describe  the  orbital  motion  about  the  earth. 

Once  the  secular  effects  were  defined  in  terms  of  the  ROE,  a  nonlinear  approach 
was  developed  to  compute  adjustments  in  the  ROE,  which  reduced  the  effects  of  the 
secular  drift. 
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II.  Literature  Review 


Although  the  approach  investigated  in  this  work  provides  new  insights  into  the 
effects  of  h  perturbations  described  by  ROE,  it  is  by  no  means  the  first  attempt  to 
understand  and  minimize  the  effect.  In  previous  literature,  the  method  is  often  referred  to 
as  establishing  a  “J2  Invariant  Formation”,  a  term  first  coined  by  Schaub  and 
Alfriend  (6:78).  Many  others  including  Wiesel,  Breger  and  How,  and  Sabitini  and 
Tragesserhave  also  added  to  the  literature  on  this  problem. 

J2  Invariance  Method 

One  method  to  solve  for  B  invariant  orbital  elements  has  been  established  by 
Hanspeter  Schaub  and  Kyle  T.  Alfriend  (6:77-85).  The  method  detennined  was  based  on 
momentum  equations  found  from  the  first-order  expansion  of  the  Hamiltonian  with 
respect  to  B-  The  system  produced  two  equations  that  are  functions  of  the  chiefs 
momentum  elements.  These  were  then  converted  back  into  classical  orbital  elements, 
resulting  in  expressions  that  were  functions  of  the  chiefs  classical  orbital  elements  and 
three  orbital  element  differences:  Sa,  Si,  and  Se.  This  allows  for  one  of  the  desired 
values  to  be  fixed  and  a  closed  form  solution  obtained  for  the  other  two. 

Though  a  solution  is  produced,  the  method  does  not  allow  for  a  significant 
amount  of  physical  interpretation  due  to  the  solution  set  existing  in  momentum  space.  In 
addition,  though  this  solution  was  based  on  a  non-circular  chief,  the  solution  allows  for  a 
secular  drift  in  the  differences  in  the  argument  of  perigee.  In  the  event  that  the  arguments 
of  perigee  become  out  of  phase  for  a  moderately  eccentric  chief  and  deputy,  it  becomes 
clear  that  the  relative  positions  could  vary  greatly  from  the  desired  conditions.  Also,  the 
solution  includes  the  tangent  of  the  chiefs  inclination.  Therefore,  a  singularity  arises  for 
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nearly  polar  orbits.  Although  a  perfectly  polar  orbit  cannot  be  obtained,  a  solution  still 
exists;  however,  this  solution  requires  considerably  large  contributions  from  to  the  other 
components. 

Partial  J2  Invariance 

Another  method  that  approached  partial  J2  invariance  was  presented  by  Breger 
and  How  (1:1-9).  The  method  refers  to  “Partial  Invariance”  due  to  considerable 
restrictions  needing  to  be  imposed  to  achieve  total  invariance.  The  used  a  state 
propagation  matrix,  a  matrix  that  is  a  function  of  the  initial  time  and  a  given  final  time. 
When  multiplying  a  state  transition  matrix  with  the  array  of  initial  conditions,  the  product 
is  the  array  of  conditions  at  the  final  time.  The  state  transition  matrix  used  by  Breger  and 
How  included  the  effects  of  J2  perturbations  in  the  form  of  osculating  elements.  The  drift 
produced  was  then  multiplied  by  a  weighted  norm,  allowing  for  a  definition  of  more  and 
less  favorable  drift  components.  The  method  also  included  a  contribution  due  to  the 
amount  of  fuel  required  to  overcome  the  secular  drift,  once  again  allowing  for  the 
multiplication  of  a  weighing  value.  Both  of  these  values  were  used  to  create  a  cost 
function.  This  value  was  then  optimized,  resulting  in  the  most  desirable  initial  conditions 
for  the  given  weighed  values. 

Special  Inclinations  from  Genetic  Algorithm 

In  addition,  the  results  of  this  J2  invariant  method  were  used  by  Sabitini,  et  al 
(5:97-100).  The  method  used  a  genetic  algorithm  with  the  refined  initial  conditions, 
established  through  use  of  the  method  determined  by  Alfriend  and  Schaub.  The 
algorithm  spanned  the  range  of  inclinations  and  used  a  value  function  to  determine  the 
method’s  accuracy  at  canceling  out  the  secular  drifts  after  100  orbits.  The  paper 
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concluded  that  the  method  was  fairly  consistent  and  was  able  to  reduce  the  secular  effects 
of  perturbation  to  approximately  two  thirds  to  one  percent  of  the  initial  formation; 
however,  for  two  angles  and  their  complements  the  method  was  significantly  more 
effective: 

/=63.4  ;  1 16.6 

.  _  o  o 

i=49  ;  131 

The  first  two  angles  are  known  as  the  critical  inclinations,  the  inclinations  where 
the  secular  drift  of  the  argument  of  perigee  is  zero.  The  second  pair  has  been  labeled  by 
the  author  as  “special  inclinations”  continuing  to  reserve  the  previously  established  title 
of  critical  inclination  for  the  previous  pair. 

Active  Control 

One  method  of  overcoming  the  secular  effects  of  perturbations  was  approached 
by  William  Weisel  (10:74-78).  The  method  included  optimizing  two-impulse  control 
maneuvers  on  satellite  fonnations  by  allowing  the  satellite  to  follow  more  natural 
dynamics.  The  method  used  a  Floquet  method  of  propagation  that  included  fourteenth 
order  geopotential  as  well  as  air  drag  and  used  the  natural  perturbation  forces  to  aid  in 
minimizing  the  amount  of  impulse  required  to  maintain  the  cluster.  However,  the  paper 
stated  that  the  fuel  required  to  maintain  specific  geometries  would  far  exceed  practicality 
and  looked  more  at  formations  intent  on  maintaining  relative  distances.  To  contrast,  this 
thesis  will  assume  that  certain  aspects  of  the  relative  geometry  are  required  and  will  try  to 
obtain  minimum  drift  rates  while  achieving  those  specifications. 
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Defining  a  Relative  Orbit 


The  method  of  obtaining  relative  geometry  between  satellites  starts  with  the 
declaration  of  the  primary  satellite.  This  satellite  will  be  labeled  the  “chief’  and  the  other 
satellites  will  be  considered  “deputies.”  For  this  thesis,  a  single  deputy  was  considered; 
however,  the  concepts  can  be  applied  for  multiple  satellite  formations.  When  describing 
the  relative  position,  the  most  commonly  used  coordinate  frame  is  Hill’s  frame.  This 
frame  is  an  orthogonal  coordinate  system  whose  origin  is  located  at  the  chiefs  position. 
The  x-axis  is  in  the  radial  direction  of  the  chief.  The  y-axis  is  in  the  direction  of  the 
chiefs  velocity,  and  the  z-axis  is  nonnal  to  the  chiefs  orbital  plane  with  a  right-handed 
orientation.  In  this  coordinate  frame  a  linear  set  of  differential  equations  describing  the 
relative  motion  was  established  by  Hill  as  well  as  Clohessy  and  Wiltshire,  with  a  slightly 
different  orientation.  Because  these  equations  are  much  alike,  the  differential  equations 
established  are  referred  to  as  Hill’s  equations,  Clohessy-Wiltshire  (CW)  equations,  and 
Hill-Clohessy-Wiltshire  (HCW)  equations.  In  the  simplest  fonn  these  differential 
equations  can  be  expressed  as: 

x-2ny -3n2x  =  0 

y  +  2nx  =  0  (1) 

z  +  n2z  =  0 

Where  n  is  the  mean  motion  of  the  chiefs  orbit  and  can  be  expressed  as: 


This  particular  form  of  these  differential  equations  is  based  on  three  assumptions. 
One  is  that  the  only  acceleration  is  gravitational  acceleration  from  a  point  mass.  The 
second  is  the  assumption  of  a  circular  chief  orbit.  The  third  is  that  the  relative  vector  is 
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considerably  smaller  than  the  chiefs  position  vector,  allowing  the  second  order  terms  to 
be  neglected  (9:282-285): 


e=0 


(3) 


It  should  be  pointed  out  that  these  differential  equations  are  independent  of  the 
deputy’s  position  in  y-direction.  This  means  that  for  unperturbed  motion,  there  are 
certain  significant  initial  conditions: 


X  = 


X 

"0" 

y 

To 

z 

0 

X 

0 

y 

0 

z 

0 

(4) 


The  resulting  relative  motion  for  these  initial  conditions  would  be  unchanged  and  has 
been  used  repeatedly  due  to  its  stable  relative  position.  This  particular  formation  is 
commonly  referred  to  as  the  leader-follower  formation.  In  addition  to  being  invariant  to 
displacements  in  the  y-direction,  the  equations  are  also  decoupled  in  the  z-direction. 
Consider  the  following  initial  conditions: 


X 

“0" 

y 

To 

z 

zo 

X 

0 

y 

0 

z 

>0. 

(5) 


This  would  produce  a  relative  orbit  that  is  invariant  in  the  x  and  y  directions,  while 
oscillating  in  the  z-direction.  Though  it  is  possible  for  large  oscillations  in  the  z-direction 
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to  make  contributions  in  the  negative  x  position,  due  to  the  deputy’s  position  projected 
onto  the  chiefs  orbital  plane,  these  linear  differential  equations  do  not  account  for  this. 

Defining  ROE 

Though  Hill’s  frame  coordinates  accurately  describe  the  relative  position  and 
velocity,  much  like  Earth  Centered  Inertial  (ECI)  position  and  velocity  vectors  for  a 
single  satellite,  these  usually  do  not  allow  for  a  very  good  understanding  of  the  relative 
orbit.  To  correct  this  inconvenience,  the  Relative  Orbital  Elements  (ROE)  were 
introduced  by  Lovell  and  Tragessor  (3:2-4)  .  These  orbital  elements  are  defined  by  the 
following  six  equations: 

ae=  2 

xd  =  4x  +  2  — 
n 

n.X 

yd=y- 2-  (6) 

n  v  ' 

P  =  atan2  ( x,  3  nx  +  2  y ) 

^rnax 

/  =  atan2  (nz,z)~  P 

Where  atan2  is  an  inverse  tangent  function  that  puts  the  first  term  over  the  second  term,  it 
also  uses  the  signs  of  each  to  eliminate  quadrant  ambiguity. 

For  a  circular  chief  and  a  deputy  with  a  matching  period,  the  deputy  would  follow 
a  two-by-one  ellipse  whose  semi-major  axis  is  in  the  direction  of  the  chiefs  velocity 
vector  and  is  designated  by,  ae.  The  center  of  the  two-by-one  ellipse  is  located  on  the 
point  {xd,  y<i).  The  angle  /?,  represents  the  deputy’s  position  on  the  ellipse.  It  is  defined  as 
the  angle  from  the  negative  x-axis  to  the  deputy’s  position  projected  on  the  chiefs  x-y 
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plane.  For  a  circular  chief  this  would  correspond  to  the  deputy’s  mean  angle  from  its 
own  perigee.  The  deputy’s  maximum  displacement  in  the  z-direction  defines  the  element 
Zmax ■  The  angle  y  is  sometimes  referred  to  as  the  relative  ascending  node,  but  is  in  fact 
opposite  in  sign  to  that  angle.  This  angle  represents  the  difference  in  the  oscillation  angle 
in  the  z-direction  and  the  angle  [>.  For  further  clarification,  please  refer  to  Figure  29  in 
Appendix  B. 


Linear  Approach 

One  approach,  to  minimize  the  secular  effects  of  J2  perturbations,  was  used  by 
Tragessor  and  Skrehart  (8:1-15)  and  is  one  that  I  will  refer  to  as  the  linear  approach.  The 
time  rate  of  change  of  the  orbital  element  differences  was  written  in  relative  orbital 
parameters  resulting  in  the  following  equations: 


£(M  +  d))  =  ~— 
2  a 


+  - 


C 


3  n  +  7C(1  -  |sin2  i)rj  +  7C(|cos2  i  - 


3(1-— sin2  z)  +  —  (— cos2  i-— ) 
2  77  2  2 


X; 


8  a' 

Q 

- cos  z  sin  z'(3/z  +  5)zmax  cos[0-(y  +  /?)] 


1C  c  c 

dtl  =  —cos  i  xd  -  - — J  cos  i  ae2  +  —  sin  i  zmax  cos[^  -(/  +  fi)\ 
2  a  2rja~  a 


(7) 


(8) 


Where  q  and  C  are  defined  as: 


n  2x1/2 

r/  =  (l-e  ) 

(9) 

3J2Re2n 

2a2r/4 

(10) 

Their  previous  work  concentrated  primarily  on  fixing  the  size  of  the  relative  orbit, 
ae,  and  solving  for  the  resulting  semi-major  axis  shift,  xj,  and  out  of  planar  motion,  zmax 
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resulting  in  both  of  the  equations  producing  zero  secular  changes.  The  resulting  solution 
was  incredibly  effective  at  establishing  invariance  for  co-planar  orbits.  However,  when 
out  of  planar  motion  was  introduced,  the  method  was  much  less  successful.  By 
attempting  to  choose  both  a  size  ae  and  fixing  a  magnitude  of  the  out  of  planar  motion, 
zmax,  the  two  equations  became  over-constrained  with  the  only  apparent  remaining 
variable  being  the  difference  in  semi -major  axis,  xj.  This  would  suggest  that  the  system 
of  equations  had  no  solution.  Essentially,  this  required  zeroing  out  two  independent 
equations  while  only  being  able  to  adjust  one  variable.  Upon  removing  the  constraint  on 
ae,  the  system  was  still  considerably  constrained  due  to  both  equations’  relative 
insensitivity  to  changes  in  ae. 

When  an  orbit  was  propagated  for  ten  orbits  without  modification  to  the  initial 
conditions,  the  corresponding  in-plane  drift  was  approximately  fourteen  meters  per  orbit. 
With  an  out  of  plane  drift  of  approximately  5  cm  per  orbit,  the  decision  was  made  to  use 
period  matching  to  cancel  the  drift  in  the  orbital  direction.  This  reduced  the  drift  in  the  y 
direction  to  1 .4  meters  per  orbit. 
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III.  Methodology 


This  work  began  by  establishing  how  the  initial  conditions  are  defined. 

Afterwards  the  secular  effects  of  the  h  perturbation  on  the  ROE  were  determined,  and  the 
conditions  for  establishing  invariance  were  defined.  Then,  the  linear  method  presented 
Tragesser  and  Skrehart  was  modified  and  analyzed.  Finally,  a  non-linear  method  was 
introduced,  and  the  accuracy  was  determined  through  numerical  integration. 

Determining  Initial  Conditions 

During  the  course  of  this  thesis  the  initial  conditions  for  the  chief  will  be  set. 
These  initial  conditions  will  be  given  in  classical  orbital  elements: 


Schaub  (7:606-61 1)6.  Schaub,  H.  and  Alfriend,  K.  T.,  “J2  Invariant  Relative  Orbits  for 
Spacecraft  Formations,”  Flight  Mechanics  Symposium,  (Goddard  Space  Flight 
Center,  Greenbelt,  Maryland),  January  18-20,  2002,  Paper  No.  1 1,  pp.  77-95 

7  currently  has  a  method  for  converting  orbital  elements  and  orbital  element  differences 

into  a  relative  position  vector  in  Hill’s  frame.  The  process  begins  by  defining  two  arrays 

that  contain  the  chief  and  deputy’s  initial  conditions  in  the  following  orbital  elements: 

a 

e 

(12) 

<h 

q2 

Q 
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Where  the  elements  qi,  q2,  and  6  have  been  used  to  cancel  out  the  singularity  in  the 
argument  of  perigee  for  a  zero  eccentricity  orbit,  and  have  been  defined  as: 

qx  =  ecos  a> 

q2=esma>  (13) 

0  =  co  +  v 

Due  to  small  eccentricities  the  true  anomaly,  n,  will  be  approximated  with  the  mean 
anomaly,  M.  The  symbol  §eoe  in  Eq.  (16)  would  represent  the  array  containing  the 
differences  in  orbital  elements  between  the  deputy  and  the  chief: 


5e 


—  V  —  V 

oe  doe  coe 


The  values  in  this  vector  are  then  described  as  follows: 


5e„„  = 


da 

59 

Si 

5q{ 

Sq2 

m 


(14) 


(15) 


This  conversion  is  in  the  form  of: 

=  (16) 
The  transfonnation  matrix  [A]  is  given  with  perturbations,  by  Gim  (2:962).  A  version 
without  the  perturbation  terms  can  be  found  in  Luck  (4: 1)  and  is  seen  in  Appendix  A. 

The  vector  output  of  this  expression  is  a  6x1  array  that  contains  both  the  relative  position 
and  velocity  in  Hill’s  frame: 
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X 


y 


x  = 


z 

X 


y 


z, 


(17) 


Though  the  state  transition  matrix,  [A],  allows  for  non-circular  chief  orbits,  the 

complexity  of  the  equations  reduces  the  understanding  of  the  resulting  equations. 

Therefore,  the  transformation  matrix  was  significantly  simplified  with  the  assumption  of 

a  circular  chief  orbit.  The  resulting  vectors  in  Hill’s  frame  are  as  follows: 

x  =  Sa-a  cos  05qx  -  a  sin  OSq2 
y  =  aSO  +  acosiSD 
z  =  a  sin  OSi  -  a  cos  6  sin  iSD 

x  =  ansin05qx  -ancosOSq2  (^) 

3  n 

y  =  — —5a  +  2an  cos  QSqx  +2an  sin  95q2 
z  =  an  cos  OSi  +  an  sin  i  sin  OSD, 


James  Luck  (4:1-4)  used  these  equations  and  substituted  them  into  Eq.  (6)  to  produce  the 
following  expressions,  which  can  be  used  to  transform  small  changes  in  COE  to  ROE: 

ae  =  2a^(Sqx)2  +(Sq2)2 
xd  =  5a 

yd  =  aSO  +  a  cos  iSD  -  2a  (sin  05 qx  -  cos  05 q2 )  ^  ^ 

zmax  =  a  j  (Si)2  +  (sin  iSD)2 
(5  =  atan2(sin  05 qx  -  cos  0Sq2 ,  cos  05q]  +sin  05 q2  ) 
y  =  atan2(sin  OSi  -  cos  0  sin  iSD,  cos  OSi  +  sin  0  sin  iSD)  -  f5 

These  expressions  may  then  be  used  to  determine  the  contributions  due  to  the 
secular  effects  of  the  orbital  elements  and  orbital  element  differences. 
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Defining  ROE  Invariance 


Though  Eq.  (19)  was  established  with  an  orbital  element  s  set  specifically 

designed  to  minimize  singularities  near  zero  eccentricity,  it  is  easier  to  understand  the 

perturbation  effects  when  the  expressions  are  written  in  classical  orbital  elements  and 

classical  orbital  element  differences.  After  simplifying  the  expressions  to  reflect  classical 

orbital  elements  and  orbital  element  differences  the  equations  become: 

ae  =  2  a8e 
xd  =  8a 

yd  =  aSO  +  a  cos  i8D.  -  2a8esin(9  -cod) 

zmax  =  a^(Si)2  +  (sin  i8Q)2  (20) 

/ 3  =  atan2  (sin(0  -  cod  ),  cos(0  -  cod  ))  =  {0  -  cod  ) 
y  =  atan2  (sin  dSi  -  cos  9  sin  fiXT,  cos  6Si  +  sin  6  sin  fiXT)  -  / 3 

It  is  clear  to  see  that  one  of  the  deputy’s  classical  orbital  elements,  the  argument 
of  perigee,  remains  in  the  equations;  fortunately,  this  allows  for  an  easier  understanding 
of  the  secular  drift  of  the  ROE  containing  that  element. 

It  is  clear  to  see  that  with  the  circular  chief  the  ROE  ae  and  xc{  will  not  experience 
secular  drift  due  to  purely  periodic  perturbations  in  the  semi-major  axis  and 
eccentricities.  In  the  equation  for  v*  the  periodic  effects  of  aS9  and  -2a8esm(0  -  cod) 
are  equal  and  opposite  for  unperturbed  dynamics.  This  is  easier  to  spot  when 
replacing aS0  +  a cos i 8D.  with y  and  by  replacing  - 2a Se sin(0 - cod)  with-ae  sin(/?) . 

yd=y~ae  sin  A  (21) 

Therefore,  in  order  to  cancel  out  secular  perturbations  in  yd  it  is  essential  to  minimize  the 
drift  rates  in  y.  It  can  be  seen  that  the  only  secularly  affected  components  of  y  are  due  to 
secular  contributions  in  the  form  of  80  and  dtl .  Therefore,  by  matching  the  secular  drift 
rates  in  the  argument  of  latitude  and  the  longitude  of  the  ascending  node,  this  ROE  will 
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not  vary  secularly.  In  addition,  it  can  be  seen  that  zmax  is  also  a  function  of  dtl  and  that 
minimizing  this  will  also  establish  reduced  secular  effect  in  zmax.  This  has  been  shown 
time  and  again  for  invariant  formations  and  could  be  expected. 

In  addition,  in  order  to  establish  true  invariance  in  the  ROE,  other  terms  would 
also  have  to  be  taken  into  consideration.  Canceling  out  secular  perturbations  found  in 
P  and  y  would  require  a  different  definition.  Due  to  the  rapid  rate  of  change  of  P ,  an 
invariant  drift  rate  could  be  defined  by  stating  a  desired  P  for  a  given  angle  0  .  In  order 
to  maintain  this  desired  initial  condition,  the  secular  drift  of  the  deputy’s  argument  of 
perigee  would  have  to  be  zero.  This  is  possible,  but  would  require  the  deputy  to  be 
located  at  the  critical  inclination. 

If  the  secular  drift  in  P  would  be  considered  zero  when  there  are  no  longer  any 
affects  of  J2,  then  the  following  condition  would  have  to  be  met: 

(M  c+d>c)  =  d>d  (22) 

It  can  be  seen  that  in  order  for  this  and  SO  to  both  be  zero,  the  secular  drift  rate  of  the 
deputy’s  mean  anomaly  would  have  to  be  zero.  This  would  also  impose  a  restriction  on 

o 

the  deputy’s  inclination;  now  requiring  a  54.7  .  Although  both  of  these  would  be 
possible,  the  imposed  inclination  restrictions  would  severely  restrict  the  orbits. 

Therefore,  canceling  out  the  secular  effects  in  the  angle  P  will  not  be  required  to 
establish  invariance. 

Furthermore,  because  y  is  equal  to  a  function  of  XI  minus  P ;  once  Stl  is  near 
zero  the  value  for  y  will  have  a  rate  of  change  equal  and  opposite  to  that  of  P .  Therefore, 
with  a  circular  chief  a  J2  invariant  relative  orbit  will  require  zero  secular  drift  i  n  yd  and 

Zmax • 
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Defining  Initial  Conditions 


The  purpose  of  this  thesis  is  to  take  a  set  of  initial  conditions  for  the  chief  and  a 
desired  set  of  ROE  and  adjust  the  ROE  slightly  in  order  to  reduce  the  secular  effects  of 
the  J2  perturbation.  The  method  will  then  be  verified  with  numerical  integration.  The 
verification  process  will  begin  by  using  the  ROE  to  determine  the  orbital  element 
differences.  The  classical  orbital  elements  for  the  deputy  are  then  calculated.  Both  sets 
of  orbital  elements  are  then  individually  transferred  into  vectors  in  the  Earth  Centered 
Inertial  (ECI)  frame.  These  were  then  integrated  using  a  fourth  order  numerical 
integrator.  The  output  was  then  converted  back  into  classical  orbital  elements,  and  the 
relative  vectors  were  then  calculated  in  Hill’s  frame.  The  ROE  were  then  calculated  for 
each  time  step. 

Due  to  the  effects  of  the  h  perturbation  being  inversely  proportional  to  the  semi¬ 
major  axis,  satellites  in  Low  Earth  Orbit  will  be  analyzed.  The  primary  semi-major  axis 
used  in  this  thesis  is  7000  km.  This  orbit  should  be  small  enough  to  produce  non-trivial 
effects  of  J2,  while  still  remaining  high  enough  that  the  effects  of  J2  will  remain  the 
dominant  perturbation.  The  first  case  will  consider  a  chief  orbit  with  an  inclination  of 
thirty  degrees.  This  inclination  has  been  chosen  because  it  is  near  the  28.5°  launch 
latitude  of  Cape  Canaveral  and  it  has  not  shown  any  traits  that  would  render  it  an 
exceptional  case.  This  case  takes  into  consideration  both  the  effects  due  to  changes  in 
eccentricity,  as  well  as  out  of  planar  motion.  The  unmodified  set  of  initial  conditions  is 
as  follows: 
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Chief  orbit 

Desired  ROE 

a  =  7000  km 

ae  = 500  m 

e  =  0 

xd  =  0  m 

i  =  30 

yd=  o  m 

[0 

ll 

o 

ll 

o 

co  =  0° 

o 

II 

v  =  0° 

Zmax=500m 

Figure  1  is  a  plot  of  the  chiefs  position,  (0,0),  and  the  deputy’s  relative  orbit  about  that 
position.  In  the  graph  the  nearly  horizontal  line  represents  the  direction  of  the  orbit 
travel.  The  nearly  vertical  line  points  in  the  negative  x  direction  representing  the 
direction  of  the  center  of  the  chiefs  orbit  and  terminates  at  approximately  the  deputy’s 
initial  position.  In  Figure  2  the  chief  is  again  located  at  (0,  0)  and  the  deputy’s  projection 
onto  the  x-y  orbital  plane  follows  the  outside  trajectory  beginning  at  (0,  -250m)  and 
proceeding  in  the  anti-clockwise  direction.  Also  plotted,  beginning  at  (0,  0)  and  drifting 
in  the  negative  y  direction,  is  xd  vs.  yd. 


Figure  1:  3-Dimensional  J2  Drift  for  Uncorrected  IC 
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Deputy  vs  Chief  and  xd  vs  yd 


Y(m) 


Figure  2:  Orbital  Plane  J2  Drift  for  Uncorrected  IC 

The  figures  above  allow  for  better  understanding  of  the  relative  motion,  as  well  as 
a  slightly  better  understanding  of  the  effects  of  the  secular  drifts  and.  In  Figure  2  the  x 
andy  components  of  the  data  have  been  plotted  in  the  chiefs  orbital  plane.  It  is  apparent 
in  this  graph  that  this  method  produces  the  classic  two-by-one  ellipse  that  defines  a 
relative  orbit  about  a  circular  chief  orbit.  Though  it  is  difficult  to  determine  the 
magnitude  of  the  secular  drift  from  this  graph,  it  is  apparent  that  a  secular  drift  is  present 
and  the  plot  of  (xd,Vd)  shows  that  the  deputy  is  drifting  in  the  negative  y-direction.  To 
allow  for  a  better  determination  of  the  rates  of  change,  the  following  figures  plot  the  y- 
position  in  meters  vs.  time  in  orbits. 
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Y  position  vs  Time 


Figure  3:  Relative  Y-Position  vs.  Time  for  Uncorrected  IC 


Y  position  vs  Time 


Figure  4:  Maximum  Displacement  in  Y-Direction  for  Uneorrected  IC 

From  Figure  3  it  is  apparent  that  the  primary  motion  in  the  y-direction  is  a 
periodic  oscillation.  However,  because  of  the  scale  it  is  hard  to  distinguish  the  secular 
effects  that  have  been  introduced  due  to  perturbation  accelerations.  Figure  4  is  a  close  up 
of  the  maximum  displacements  in  the  v-direction.  The  figure  shows  that  there  is  a 


20 


negative  drift  of  approximately  negative  one  meter  per  orbit.  Likewise  the  following  two 
graphs  plot  the  displacement  in  the  z-direction  vs.  time. 


Z  Position  vs  Time 
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Figure  5:  Out  of  Plane  Position  vs.  Time 


Z  Position  vs  Time 


Figure  6:  Maximum  out  of  Plane  Position  vs.  Time 

Once  again,  Figure  5  shows  a  dominant  harmonic  oscillation  in  the  z-direction; 
however,  the  scale  in  Figure  6  allows  for  a  better  understanding  of  the  drift  in  the  z- 
direction.  This  drift  is  nonlinear,  increasing  about  70  cm  over  25  orbits. 
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Pure  Out-of-Plane  Relative  Motion 


In  the  event  that  the  relative  motion  desired  is  a  pure  oscillation  in  the  z-direction, 
the  value  of  ae  will  be  set  to  zero.  This  will  cause  both  /?  and  y  to  lose  definition,  as  both 
become  a  function  of  zero  divided  by  zero.  However,  because  this  singularity  exists  in 
the  transformation  from  Hill’s  Frame  components  into  ROE,  declaration  of  these  values 
will  produce  legitimate  initial  conditions  in  Hill’s  frame  that  can  be  propagated  to 
determine  the  relative  motion. 

This  relative  motion  would  include  a  hannonic  oscillation  in  the  z-direction  that 
intersects  the  chiefs  orbital  plane  at  (yd,  Xd).  Therefore,  in  order  for  the  initial  conditions 
to  be  more  realistic  this  case  will  include  a  non-zero  Vd  to  prevent  a  collision  with  the 
chief  satellite.  However,  due  to  previously  established  independence  to  displacements  in 
the  y-direction  this  will  not  affect  the  stability  of  the  relative  orbit. 

The  initial  conditions  have  been  established  as  follows: 


Chief  orbit 

Desired  ROE 

a  =  7000  km 

at=0m 

e  =  0 

xd  =  0  m 

i  =  30° 

yd  = -25  m 

Q  =  0° 

ll 

o 

<x>  =  0° 

y  =  0° 

v  =  0° 

Zmax  =  500  m 

These  initial  conditions  were  numerically  propagated  and  can  be  seen  in  Figures  7 

and  8. 
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Y  Position  vs  Time 
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Figure  7:  Maximum  Displacement  in  Y-Direction 


Z  Position  vs  Time 


Figure  8:  Maximum  out  of  Plane  Position  vs.  Time 

In  Figure  7,  it  can  be  seen  that  it  is  no  longer  necessary  to  scale  the  figure  in  order 
to  determine  the  secular  drift  in  the  v-dircction,  due  to  the  absence  of  the  large  harmonic 
oscillations  that  were  present  for  ae= 500  m.  This  figure  also  shows  a  significant  amount 
of  secular  drift  in  the  orbital  direction,  approximately  -7  meters  per  orbit. 
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Figures  6  and  8  show  that  in  both  cases,  the  drift  in  the  z-direction  appears 
parabolic.  This  is  partly  because  the  primary  contribution  to  zmax  is  due  to  a  difference  in 
inclinations.  Given  the  time  derivative  of  zmax: 


d (zmax )  _  a  sin2  iSClStl 
dt  ^(Sif  +  (sin  iXT)1 

For  ( Si)2»  (sin /TO),  the  equation  is  approximately  equal  to: 

d(z max)  <2 sin2  iSClStl 

dt  \Si\ 


(23) 


(24) 


With  a  relatively  constant  value  of  the  value  SO-  would  increase  linearly. 
Therefore  Eq.  (24)  would  increase  linearly,  resulting  in  parabolic  growth  of  zmax.  This 
would  be  indicative  of  a  chief  orbit  with  fairly  low  inclination,  in  which  much  of  the  out 
of  plane  motion  is  due  to  the  difference  in  inclination  as  opposed  to  the  right  ascension  of 
the  ascending  node.  This  is  the  case  with  the  two  previous  scenarios,  and  helps  explain 
the  growth  of  zmca  in  Figures  6  and  8. 

In  the  event  that  (sin  iSOl)2»  (Si)2  the  equation  can  be  approximated  as: 


^Zmax^  «  a  sin  i Stl 
dt 


(25) 


This  equation  would  produce  linear  drift  that  would  be  expected  due  to  secular  effects. 
This  would  represent  a  chief  orbit  that  is  near  polar,  where  most  of  the  out-of-plane 
motion  is  due  to  difference  in  the  right  ascension  of  the  ascending  node.  In  both  cases,  a 
zero  secular  drift  in  the  difference  in  the  longitude  of  the  ascending  node  would  result  in 
zero  rate  of  change  of  zmax.  Therefore,  the  secular  drift  rate  of  the  difference  in  the 
longitude  of  the  ascending  nodes  will  be  analyzed  in  more  depth. 
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Secular  Expansion 


The  equation  for  the  secular  drift  in  the  differences  of  the  longitude  of  the 
ascending  node  is  as  follows: 


m  = 


3  J2Rend 


2“d  (l~ed  ) 


2,2  C0SZ^ 


2J2R;nc 

2ac2(l-ec2)2 


cosz.. 


(26) 


After  including  the  assumption  of  circular  chief  orbit  and  writing  strictly  as  a  function  of 
the  chiefs  orbital  elements  and  the  orbital  element  differences  the  equation  becomes: 


m  = 


u2Re2M 1/2 


2{a  +  Sa)1 12  Se1)2 


cos  (i  +  di) 


'  '  3 J,R,2m"2 


2  a 


7/2 


cosz 


(27) 


Because  the  orbital  element  differences  are  considerably  smaller  than  the  terms  they  are 
grouped  with,  the  expressions  containing  the  differences  can  be  expanded  about  the 
dominant  value.  After  performing  the  expansions  and  combining  the  higher  order  terms, 
the  expression  becomes: 
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2 a 


7/2 


In  order  to  better  understand  the  accuracy  of  the  expression  it  would  be 
worthwhile  to  consider  the  approximate  order  of  magnitude  of  the  higher  order  terms. 
For  a  semi-major  axis  of  around  7000  km,  in  order  to  produce  a  zmax  equal  to  or  less  than 
5km  the  maximum  Si  will  be  on  the  order  of  10‘4.  Likewise,  for  an  ae  on  the  order  of 
lkm,  the  resulting  Se  will  be  on  the  order  of  10'\  Also,  because  secular  drift  due  to 
differences  in  mean  motion  is  highly  dependent  on  differences  in  semi-major  axis,  J2 
invariance  will  often  require  Sa  smaller  than  10'4  km.  Therefore,  higher  order  terms 
inside  the  parenthesis  will  be  on  the  order  of  10'  “  radians  per  second.  The  resulting 
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round  off  error  would  create  a  secular  drift  of  zmax  on  the  order  of  millimeters  per  orbit. 
After  cancellation  and  prioritizing  by  approximate  maximum  orders  of  magnitude  the 
equation  becomes: 


Once  again,  if  a  1  km  zmax  was  created  strictly  with  a  change  in  inclination,  the 
first  tenn  could  be  on  the  order  of  10‘4,  whereas  the  remaining  terms  would  be 
approximately  four  orders  of  magnitude  smaller.  This  suggests  that  the  largest  step 
towards  establishing  E  invariance  of  this  orbital  element  difference  is  to  minimize  the 
difference  in  the  inclinations.  However,  Eq.  (20)  shows  that  the  magnitude  of  the  out  of 
planar  motion,  zmax,  is  created  by  a  combination  of  Si  and  SQsini.  Therefore,  for  any 
chief  orbit  without  a  near  zero  inclination  this  would  not  limit  zmax ;  this  would  only  limit 
the  contribution  resulting  from  the  difference  in  inclinations. 


Linear  Approach 

With  this  idea  in  mind,  the  linearizations  put  forth  by  Tragesser  and  Skrehart  were  then 
taken  into  consideration.  Previous  work  has  shown  that  these  produce  practical  results 
for  co-planar  orbits;  however,  when  out-of-planar  motion  is  introduced  the  method 
becomes  much  less  successful  due  to  the  sensitivities  in  the  equation  to  small  changes  in 
zmax.  However,  if  it  was  possible  to  nearly  zero  out  the  coefficients  of  zmax  the  equations 
would  become  much  less  sensitive  to  out-of-plane  motion  and  it  would  be  possible  to  re¬ 
establish  the  stability  shown  for  co-planar  orbits.  Therefore,  the  coefficients  of  zmax  will 
be  considered  in  more  depth: 
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(30) 


C 

S(M  +  d>)  =  fl(Coe,xd,ae  ) - cos/sin  i(3n  +  5 )zmax  cos  [9  -  (y  +  f3 )] 

a 

SCI  =  f2(Coe,xd ,  ae2 )  +  —sin  i  zmax  cos[9  -  (y  +  /?)] 

a 

Two  terms  that  are  present  in  both  coefficients  of  zmax  are  aA  and  sin/.  Though  it 
would  be  possible  to  increase  the  magnitude  of  the  variable  a  to  the  point  where  both 
equations  become  nearly  independent  of  zmax,  and  even  more  possible  to  choose  a  chief 
with  a  zero  inclination,  the  resulting  orbits  would  be  considerably  constrained.  There  is 
however  another  term  that  shows  up  in  both  equations: 

cos  [9  -(/  +  /?)]  (31) 

This  term  is  bound  between  one  and  negative  one.  In  theory,  if  it  is  possible  to 
choose  the  correct  combination  of  9,  y,  and  [>  it  should  be  possible  to  nearly  zero  out  the 
contributions  of  zmax  and  re-establish  the  same  stability  displayed  for  co-planar 
formations. 

Because  y  is  the  angular  difference  between  the  oscillation  angle  in  the  z-direction 
and  the  position-based  angle  /?,  (y  +  J3)  simplifies  to  the  angle  of  the  oscillation  in  the  z- 
direction  and  will  be  labeled  as  yz. 

yz  =  atan2(«z,z)  =  y  +  fi  (32) 

cos[9-(y  +  /?)]  =  cos(9-y_)  (33) 

These  angles  left  inside  the  parenthesis  are  increasing  at  a  given  frequency.  For  yd  =  0 

both  angles  would  increase  at  a  rate  equal  to  the  mean  motion.  Resulting  difference  is 
relatively  constant.  Because  there  are  now  two  components  that  are  functions  of  relative 
orbital  elements,  it  is  possible  to  fix  znmx  and  replace  the  expression  cos[<9  -  (y  +  /?)]  with 
a  new  variable  zy. 
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zy=cos(9-yz)  (34) 

Furthermore,  if  a2  is  also  chosen,  two  of  the  main  components  of  the  relative 


geometry  have  been  fixed;  the  magnitude  of  the  two-by-one  ellipse,  a2,  and  the 
magnitude  of  the  out  of  plane  motion,  zmax.  By  setting  the  difference  in  drift  rates  equal 
to  zero  and  bringing  the  now  constant  contribution  of  ae2  to  one  side  the  expression  can 
be  written  as: 


This  gives  a  closed  form  solution  for  both  zy  and  x(i  in  the  form: 

Xd]  =  UV  [D]  ae2  (37) 

z 

V  y  J 

Note  that  since  z},  is  a  function  of  the  sum ;/  +  /?,  we  can  arbitrarily  select  /?  and 
still  satisfy  the  conditions  by  choosing  the  correct  y  to  satisfy  Eq.  (32).  This  means  that 
for  a  given  zmax  and  ae  the  full  solution  is  now  a  solution  set  whose  relative  positions  lay 
on  the  two-by-one  ellipse  offset  by  the  chiefs  orbital  plane  by  a  distance  ofzmax  sin  y, . 
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This  would  allow  for  multiple  satellites  to  occupy  the  same  solution  set,  with  different 
angles  fi  and  y  but  still  maintaining  the  same  zr  In  addition,  because  the  cosine  function 
is  even,  both  a  positive  and  negative  angle  y2  fulfill  the  requirements  for  zr  Combine  this 
with  yd  independence,  and  this  allows  for  fi,  zmax,  ae,  and  Vd  to  be  set  arbitrarily.  This 
allows  the  ROE  Xd  and  y  to  be  adjusted  in  order  to  establish  the  desired  invariance. 

Analysis  of  Linear  Approach  for  Circular  Deputy 

To  investigate  the  effectiveness  of  this  method  to  reduce  the  effects  of  the  J2 
perturbations  on  a  formation,  numerical  simulations  were  performed.  This  method  will 
include  the  initial  conditions  established  previously  for  a  pure  oscillation  in  the  z- 
direction.  Once  again,  this  initial  condition  causes  fi  and  y  to  lose  physical  interpretation 
and  become  undefined;  however,  because  zy  remains  defined,  the  initial  value  of  fi  will 
once  again  be  set  to  zero  and  the  value  for  the  initial  condition  of  y  will  be  computed. 


Chief  orbit: 


Desired  ROE: 


a  =  7000  km 
e  =  0 
i  =  30° 

0  =  0° 
co  =  0° 
v  =  0° 


fle=0m 

xd  =  0  m  +  TBD  Adjustment 
y d  = -25  m 


p  =  o° 
y  =  TBD 


After  substituting  these  values  into  the  previously  established  equations  the  initial 
conditions  are  determined  to  be: 


(x  ) 

Ad 

0  km 

\zy  ( 

0 

(38) 
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The  zy  allowed  for  a  y=  ±  90°.  For  this  case  the  positive  value  of  gamma  will  be  chosen. 
The  resulting  initial  conditions  are  as  follows: 


Chief  orbit: 

a  =  7000  Ion 
e  =  0 
i  =  30° 

Q  =  0° 
co  =  0° 
v  =  0° 


Desired  ROE 

ae  =  0  m 
Xj  =  0  mm 
yd  =  -25m 
/?  =  0° 
y  =  90° 
zmax  =  500  m 
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Figure  9:  Maximum  Displacement  in  F-direction 
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Figure  10:  Maximum  Displacement  in  Z-direction 
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As  Figure  9  shows,  this  method  is  highly  successful  at  reducing  the  secular  drift 
in  the  v-direction.  The  resulting  secular  drift  is  approximately  -1  mm/orbit  as  opposed  to 
over  -7  m/orbit  drift  seen  in  Figure  7.  Likewise,  Figure  10  also  shows  a  drift  multiple 
times  smaller  than  the  drift  observed  in  Figure  8. 

It  is  also  worth  mentioning  that  the  apparently  trivial  solution  seen  in  Eq.  (38) 
allows  for  a  slight  amount  of  understanding  of  the  modifications  to  the  initial  linear 
method.  Had  the  system  of  equations  been  used  as  they  were  previously  in  Eqs.  (7)  and 
(8),  to  solve  zmax  and  xj,  the  solution  set  would  have  required  a  zero  maximum 
displacement  in  the  z-direction.  This  would  limit  the  solution  to  a  leader-follower 
formation  observed  in  Eq.  (4).  However,  by  solving  for  zy  the  solution  now  allows  for 
displacements  in  the  z-direction,  resulting  in  the  initial  conditions  previously  mentioned 
in  Eq.  (5). 


Analysis  of  Linear  Approach  for  Non-Circular  Deputy 


Though  this  method  produces  accurate  results  for  near  zero  ae,  for  moderate  ae 
this  method  produces  much  less  desirable  results.  The  following  initial  conditions  are 
identical  to  the  previously  chosen  initial  conditions  with  the  exception  of  a  500  meter  ae : 
For  a  given  chief  orbit:  Desired  ROE: 


a  =  7000  km 
e  =  0 
i  =  30° 

0  =  0° 
co  =  0° 
v  =  0° 


ae  =  500  m 

xd  =  0  +  TBD  Correction 

yd  = 0'» 


/?  =  0° 
y  =  TBD 
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For  this  particular  formation  Eq.  (37)  establishes  xd  =  -0.03 1  mm  and  y  =  89.996° .  The 

graph  below  shows  the  projection  of  the  deputy’s  orbit  projected  onto  the  chiefs  orbital 
plane. 
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Figure  11:  Relative  Orbit  from  Linear  Results  ae  of  500  m 
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Figure  12:  Maximum  Positive  Displacement  in  F-direction  from  Linear  Results 


The  Figures  1 1  and  12  clearly  show  that  yd  is  drifting  in  the  positive  v-direction  at  a  rate 
of  over  5  meters  per  orbit.  A  second  run  with  x(j=  0m  gave  almost  identical  results. 
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When  also  considering  that  positive  displacements  in  Xd  result  in  negative  drifts,  it  can 
safely  be  deducted  that  %d  should  have  been  larger  in  magnitude  and  positive.  However, 
this  case  still  shows  that  this  method  is  still  highly  successful  at  reducing  the  secular  drift 
in  zmax,  as  shown  in  Figure  13. 
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Figure  13:  Maximum  Displacement  in  Z-Direction 
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Due  to  the  less  than  desirable  accuracy  of  the  linear  solver,  a  more  accurate  method  of 
determining  Xd  was  required. 


Single  Variable  Nonlinear  Solver 

Since  the  linear  approach  was  unable  to  cancel  out  the  secular  drift  in  the  orbital 
direction,  a  nonlinear  approach  was  taken  to  find  values  of  xd  and  y  which  would  result 
in: 


yd  = 0 

z„„„  =  0 


fiy , 

In  order  to  determine  the  initial  conditions  from  the  linearization  were 

Sx, 


numerically  propagated  for  four  orbits.  The  time  and  position  of  the  maximum 
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displacement  in  the  v-direction  were  recorded  and  run  through  a  linear  best  fit  to  help 
reduce  error.  The  slope  that  was  produced  through  best  fit  was  considered  the  secular 
drift  in  the  v-direction,  yd .  This  was  repeated  for  a  second  xd.  The  difference  in  the 


c *• 

slopes  produced  by  the  linear  best  fit,  over  the  differences  in  xd  was  and  — was 

Sx, 


determined.  The  secular  drift  yd0  was  the  error  and  the  new  xd  was  calculated: 


er  =yd  o 


yd 


—  6 


syd 

Sx, 


Sx, 


Setting  yd  equal  to  zero  and  solving  for  xd  produces: 


(39) 


5x.,=-eSx‘ 


Syd 

xd=8xd+xd  o 


(40) 


Due  to  non-linearity  in  the  equation  this  method  had  to  be  repeated  until  the  drift  rate  was 
within  desired  tolerances.  The  xd  was  then  re-introduced  into  the  secular  drift  rate  of  Stl 
to  solve  fory . 


Thirty  Degree  Inclination 


After  using  the  non-linear  approach  to  the  previous  set  of  initial  conditions 
resulted  in  the  following  relative  orbit: 

For  a  given  chief  orbit:  Desired  ROE: 


a  =  7000  km 
e  =  0 
i  =  30° 

0  =  0° 
co  =  0° 
v  =  0° 


ae  =  500  m 
xd  =  67.61  cm 

yd  =  0  rn 

/3  =  0° 
y  =  89.972° 


=  500  m 
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The  resulting  relative  orbit  is  shown  in  Figure  14. 


Deputy  vs  Chief  and  xd  vs  yd 


Y(m) 

Figure  14:  Corrected  Relative  Orbit  with  30°  Inclination  Chief 


Y  Displacement  vs  Time 


Figure  15:  Maximum  Y  Displacement  vs.  Time  with  30°  Inclination  Chief 
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Z  Position  vs  Time 


Figure  16:  Maximum  Z  Displacement  vs.  Time  with  30°  Inclination  Chief 

It  can  be  seen  from  Figures  14  and  15  that  this  method  is  successful  at 
significantly  reducing  the  secular  drift  rate  in  the  orbital  direction,  while  Figure  16  shows 
that  the  method  has  also  maintained  accuracy  in  the  out-of-plane  direction.  This  case  was 
for  verification  purposes  only;  therefore,  the  iterations  were  stopped  at  an  accuracy  of 
approximately  1  cm  drift  per  orbit  in  the  v-dircction.  When  the  method  is  used  to 
maximize  accuracy  the  iterations  will  continue  until  the  secular  drift  in  the  orbital 
direction  is  less  than  1  mm  per  orbit. 

Two  Variable  Nonlinear  Method 

Although  the  single  variable  method  is  effective  at  reducing  the  secular  effects  of 
the  h  perturbation,  it  is  possible  to  further  reduce  the  secular  drift  in  the  out-of-plane 
direction  by  expanding  the  non-linear  method  to  solve  for  both  variables  simultaneously. 
The  method  begins  with  a  propagation  of  the  initial  conditions  for  the  previously 
established  four  orbits.  The  secular  drifts  in  they  and  z  directions  are  detennined. 
Afterwards  the  obit  was  propagated  once  again  with  only  a  small  change  in  Xd  and  then 
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again  with  only  a  small  change  in  y.  The  resulting  rates  of  change  of  the  maximum 
displacements  in  the  y  and  z  directions  were  once  again  determined.  These  values  were 
then  used  to  express  the  differences  in  the  rates  of  change  of  the  error  over  the  differences 
in  the  initial  conditions. 


The  error  was  then  defined  as: 


4 vd 


A xd 


4 h 
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A^ma 
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( yd  o  ^ 


V^maxO  J 


(41) 


(42) 


The  initial  conditions  for  the  next  iteration  were  then  detennined  by  the  following 
equation: 
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(43) 


After  applying  this  method  to  the  initial  conditions  previously  established  for  the 
thirty  degree  chief  inclination,  the  following  initial  conditions  describe  the  chiefs  orbit 
and  the  relative  orbit. 

For  a  given  chief  orbit:  Desired  ROE: 

a  =  7000  km  ae  =  500  m 

e  =  0  x^=68.35cm 


i  =  30° 
Q  =  0° 
co  =  0° 
v  =  0° 


yd  =  0  rn 

4  =  0° 


y  =  90.4442° 


=  500  m 
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These  initial  conditions  were  propagated  for  twenty-five  orbits  and  the  drift  rates  in  the  y 
and  z  directions  are  shown  in  Figures  17  and  18. 


Y  Displacement  vs  Time 


Orbit  Number 

Figure  17:  Secular  Drift  in  Y-Direction  from  Two-Variable  Non-Linear  Method 


Z  Position  vs  Time 


Figure  18:  Secular  Drift  in  Z-Direction  from  Two-Variable  Non-Linear  Method 
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Figures  17  shows  that  the  secular  drift  in  the  y-direction  has  been  successfully  reduced  to 
under  5  cm  for  the  propagation  period  shown.  Figure  18  shows  that  this  method  has  been 
able  to  reduce  the  secular  drift  in  the  z-direction  to  approximately  3  cm  for  the 
propagation  period.  This  is  a  significant  improvement  over  the  approximately  20  cm  seen 
with  the  single  variable  method. 
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V.  Results  and  Discussion 


Circular  Chief 

After  verifying  that  the  single-variable  non-linear  method  was  now  successful,  a 

series  of  inclinations  were  used  to  further  analyze  this  method.  Additional  inclinations 

will  include  the  inclination  of  the  international  space  station,  the  critical  inclination,  a 

true  polar  orbit,  and  a  sun-synchronous  orbit.  The  initial  conditions  are 

Chief  Orbital  Elements:  ROE:  Chief  Inclinations: 

a  =  7000  km  ae=500m 

e  =  °  yd  =  0m 

Q  =  °;  Zmax  =  500m 

CO-0  „ 

v  =  0° 

After  applying  this  method  for  these  chief  inclinations  the  resulting  initial  conditions 
were  detennined  and  can  be  seen  in  Table  1. 


30° 

51.6° 

63.4° 

90° 

95.4° 


Inclination 

xd 

Y 

30°  Inclination 

67.2  cm 

89.582° 

International  Space  Station  (5 1 .6°) 

67.7  cm 

90.213° 

Critical  Inclination  (63.4°) 

68.1  cm 

90.173° 

Polar  Orbit  (90°) 

67.8  cm 

90° 

Sun  Synchronous  (95.4°) 

67.8  cm 

90.463° 

Table  1:  Determined  Initial  Conditions  for  Propagated  Orbits 


After  propagating  the  initial  conditions  for  the  given  cases  the  drift  in  the  z  and  y 
directions  was  plotted  in  Figures  19  and  20.  The  magnitudes  of  the  deviation  were  taken 
and  all  cases  were  plotted  together  in  order  to  allow  for  a  better  visual  comparison. 
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Magnitude  of  the  Deviation  in  the  Maximum  Value  of  Z 
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Figure  20:  Magnitude  of  the  Deviation  in  the  Maximum  Z-Direction 

Figure  19  shows  that  this  method  is  quite  successful  at  canceling  out  secular  drift  in  the 
orbital  direction,  with  the  worst  case  producing  less  than  5  cm  drift  over  the  twenty-five 
orbit  propagation  period.  Also,  with  approximately  three  centimeters  separating  the  best 
from  the  worst  case,  this  method  appears  highly  consistent  across  the  range  of 
inclinations  tested.  Figure  20  shows  that  the  ability  to  minimize  the  secular  drift  in  the  z- 
direction  is  much  less  consistent  across  the  inclinations  tested.  This  figure  does  however 
show  that  this  method  appears  highly  successful  at  inclinations  at  and  near  90°.  In  order 
to  confidently  establish  the  approximate  rates  of  change,  these  values  were  put  through  a 
linear  best  fit.  The  slopes  were  determined  in  units  of  centimeters  per  orbit  and  have 
been  listed  in  Table  2. 
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Inclination 

Drift  Rate  of  zmax  (cm/orbit) 

Drift  Rate  ofyd  (cm/orbit) 

30°  Inclination 

-1.48 

-0.13 

International  Space  Station 

0.43 

-0.19 

Critical  Inclination 

1.08 

-8.80E-02 

Polar  Orbit 

-1.02E-02 

-7.04E-02 

Sun  Synchronous 

4.11E-03 

0.12 

Table  2:  Drift  Rates  of  Maximum  Displacements  in  Z  and  Y  Directions  in  cm/orbit 


Non-Circular  Chief 


This  method  allows  for  a  J2  invariant  solution  for  two  satellites  by  establishing 
e>(M  +  of)  and  c?2both  equal  to  zero.  However,  if  it  is  possible  to  also  set  Sd>  to  zero 
then  it  might  be  possible  to  lift  the  circular  chief  assumption.  Fortunately,  if  both 
satellites  were  operating  at  the  critical  inclination  then  this  condition  would  be  met.  The 
eccentricity  of  the  chief  was  chosen  such  that  the  difference  between  the  radius  at  perigee 
and  apogee  would  be  1000  km.  Also,  the  semi-major  axis  has  been  increased  to  7500  km 
to  keep  a  safe  altitude  while  at  the  radius  of  perigee.  The  eccentricities  have  been  chosen 
to  produce  a  difference  in  the  radii  of  apogee  and  perigee  of  approximately  1000  km  and 
100  km  respectively. 

The  following  initial  conditions  have  been  chosen: 

Chief  Orbital  Elements:  Desired  ROE: 

a  =  7500  km 


i  =  63.4° 
Q  =  0° 
o-0° 
v  =  0° 


ae  =500m 

yd  =0m 


=  500m 


/?  =  0° 


The  corresponding  calculated  values  for  the  initial  conditions  are  shown  in  Table  3. 


Eccentricity 

xd 

Y 

Case  1 

0.0667 

77  cm 

90.241° 
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|  Case  2  |  0.00667  |  60.8  cm  |  90.190° 

Table  3:  Initial  Conditions  for  Non-Circular  Chief 
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Figure  21:  Magnitude  of  Drift  in  Y-Direction  for  Non-Circular  Chief  Orbits 


Magnitude  of  the  Deviation  in  the  Maximum  Value  of  Z 
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Figure  22:  Magnitude  of  Drift  in  Z-Direction  for  Non-Circular  Chief  Orbits 
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Figure  21  shows  that  the  case  of  a  moderately  eccentric  chief  is  not  as  successful 
at  canceling  the  drift  in  the  orbital  direction;  while  the  case  for  small  eccentricity  chief 
produced  results  as  accurate  as  those  produced  for  a  circular  chief.  Unfortunately,  Figure 
22  shows  that  this  method  is  nowhere  near  as  successful  at  eliminating  the  drift  in  the  z- 
direction  for  the  moderate  eccentricity  chief  orbit. 


Although  it  is  difficult  to  see  due  to  the  scale  of  Figure  22  the  drift  rate  in  the  z- 
direction  is  approximately  thrice  the  drift  seen  for  the  circular  chief  with  the  same 
inclination.  It  is  true  that  the  difference  in  semi-major  axis  is  approximately  500  km; 
however,  the  case  was  rerun  at  a  matching  semi-major  axis  and  produced  similar  results. 
In  order  to  determine  if  these  results  were  an  isolated  case,  the  same  eccentricity  of  .0067 
was  also  simulated  at  the  ISS  inclination.  The  resulting  ROE  are  given  in  Table  4. 


Eccentricity 

Inclination 

xd 

Y 

Case  1 

0 

63.4° 

68.1  cm 

90.173° 

Case  2 

0.0067 

63.4° 

68.1  cm 

90.2128° 

Case  3 

0 

51.6 

67.7  cm 

90.213° 

Case  4 

0.0067 

51.6 

70.0  cm 

90.1781° 

Table  4:  Initial  Conditions  for  Circular  and  Near  Circular  Comparisons 


The  initial  conditions  were  propagated  once  again  for  25  orbits,  producing  the  results 
shown  in  Figure  23. 
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Magnitude  of  the  Deviation  in  the  Maximum  Value  of  Z 
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Figure  23:  Deviation  in  the  Z-Direction  for  Circular  and  Near  Circular  Chief  Orbits 

It  can  be  seen  from  Figure  23  that  the  method  appears  to  produce  larger  out-of-plane  drift 
for  a  near  circular  chief  multiple  than  it  does  for  a  circular  chief. 

Additional  Drift  Rates 

With  the  previous  case  for  non-circular  chief  at  ISS  inclination,  the  initial  primary 
concern  was  the  drift  secular  drift  in  the  argument  of  perigee.  In  order  to  determine  the 
effects  due  to  secular  drift  of  the  argument  of  perigee,  the  maximum  displacement  in  the 
x-directions  were  also  compared.  It  turns  out  that  the  previous  notion  of  occupying  the 
critical  inclination  in  order  to  reduce  the  secular  drifts  in  the  difference  of  the  argument 
of  perigee  was  somewhat  effective;  however,  even  when  the  chief  occupied  the  ISS 
inclination  the  secular  drift  in  the  x-direction  was  considerably  smaller  than  the  drift  in 
the  z-direction,  as  seen  in  Figure  24. 
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Magnitude  of  the  Deviation  in  the  Maximum  Value  X 
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Figure  24:  Magnitude  of  the  Deviation  in  the  Maximum  Value  of  X  for  Non-Circular  Chiefs 

It  turns  out  that  the  drift  rate  produced  in  the  x-direction  is  not  only  small,  but  it  is 
actually  less  than  the  drift  in  the  v-dircction.  This  is  partially  due  to  the  initial  conditions 
that  established  the  difference  in  the  argument  of  perigee,  §co,  to  be  near  zero.  Applying 
this  to  a  small  angle  approximation  would  result  in  the  drift  in  the  x-direction  to  be 
approximately  on  the  order  of  the  dof.  In  an  attempt  to  maximize  the  drift  rate  in  the  x- 
direction  the  case  was  run  again  with  the  arguments  of  perigee  with  a  90  separation.  In 
this  case  the  secular  drift  for  the  25  orbit  propagation  was  approximately  7.5  cm. 
However,  the  size  of  the  relative  orbit  was  increase  drastically.  If  taking  into 
consideration  that  the  new  ae  had  increased  to  over  100  km,  it  is  easy  to  see  that  the 
percent  difference  between  the  secular  drift  and  the  relative  major  axis  is  considerably 
small.  Unfortunately,  with  such  an  increase  in  the  size  of  the  relative  orbit,  it  was  also 
concluded  that  this  method  is  considerably  less  flexibility  in  the  initial  value  of  [1. 
Therefore,  in  order  to  maintain  the  flexibility  in  the  ROE  the  case  if  a  non-circular  chief 
will  be  deemed  impractical  and  additional  work  in  this  thesis  will  only  consider  the  case 
of  the  circular  chief  reference  orbit. 
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Two-Variable  Non-Linear  Method 


After  determining  from  Figure  20  that  the  single-variable  method  was  unable  to 
produce  accurate  results  for  all  inclinations,  the  two-variable  method  was  used  at  30° 
chief  inclination.  The  results  are  plotted  in  Figure  25. 


Secular  Drift  in  Z  Direction 


Orbit  Number 

Figure  25:  Comparison  of  Secular  Drifts  in  Z  Direction 

It  is  easy  to  see  in  Figure  25  that  for  a  chief  inclination  of  30°,  the  two-variable  non-linear 
method  is  much  more  successful  at  reducing  the  secular  drift  in  the  z-direction  than  the 
single-variable  method.  Therefore,  the  method  will  be  applied  to  the  other  inclinations 
whose  secular  drifts  from  the  single-variable  method  were  still  large.  The  following 
initial  conditions  will  be  used  with  chief  orbital  inclinations  30°,  ISS  inclination,  and 
Critical  Inclination.  The  resulting  ROE  were  calculated  and  are  displayed  in  Table  5. 
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Chief  Orbital  Elements: 


Desired  ROE: 


a  =  7500  km 


a  =500m 

e  =  0 

Q  =  0° 

£ 

o 

ll 

co  =  0 

Zmax=500m 

0° 


/?  =  0° 


Inclination 

xd 

Y 

30°  Inclination 

68.35  cm 

90.444° 

International  Space  Station  (51.6°) 

67.86  cm 

90.031° 

Critical  Inclination  (63.4°) 

67.87  cm 

90.044° 

Table  5:  Initial  Conditions  for  Two-Variable  Non-Linear  Method 


After  propagating  these  initial  conditions  the  secular  drift  in  the  z-direction  was 
plotted  in  Figure  26. 


Two  Variable  Secular  Drift  in  Z-Direction 
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Figure  26:  Secular  Drift  in  Z-Direction  for  Multiple  Chief  Inclinations 

Figure  26  shows  that  this  method  is  successful  at  reducing  the  secular  drift  in  the  z- 
direction  to  fewer  than  5  cm  over  the  propagation  period  of  25  orbits.  The  secular  drift  in 
the  ^-direction  is  plotted  in  Figure  27. 
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Two  Variable  Secular  Drift  in  Z-Direction 
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Figure  27:  Secular  Drift  in  F-Direction  for  Multiple  Chief  Inclinations 

It  is  apparent  from  Figures  26  and  27  that  the  method  has  successfully  reduced  the 
secular  drifts  in  both  directions  to  the  same  order. 

After  verifying  that  the  method  has  successfully  reduced  the  secular  drift  in  the  z- 
direction  the  two-variable  method’s  initial  conditions  for  these  three  inclinations  were 
compared  with  the  single-variable  near-polar  results. 


Two  Variable  Secular  Drift  in  Z-Direction 
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Figure  28:  Comparison  of  Single-Variable  and  Two-Variable  Non-Linear  Methods 

Figure  25  shows  that  the  two-variable  method  has  produced  significantly  more 
accurate  results  than  the  single-variable  method  for  a  30°  chief  inclination.  However, 
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Figure  28  shows  that  the  method  is  still  dependant  on  the  chiefs  inclination;  with  the 
critical-inclination  being  the  only  inclination  in  this  set  that  is  able  to  produce  results  as 
accurate  as  the  near-polar  results  for  the  single-variable  method. 
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VI.  Conclusion 


This  thesis  approached  the  problem  of  establishing  J2  invariant  formations  by 
setting  the  chiefs  COE  and  stating  a  desired  set  of  ROE.  It  was  then  determined  that  the 
two  ROE,  Xd  and  y,  could  be  adjusted  in  order  to  reduce  the  secular  effects  of  this 
perturbation  acceleration.  This  work  consisted  of  establishing  a  non-linear  method  for 
determining  the  ROE  xj.  Once  this  method  had  produced  desirable  tolerances  in  the 
secular  drift  rate  in  the  orbital  direction  the  value  of  Xd  was  used  in  a  linear  function  to 
solve  for  the  angle  of  the  oscillation  in  the  z-direction.  With  this  value  and  an  arbitrary 
selection  of  the  angle  /?  the  angle  y  was  calculated,  completing  the  set  of  initial 
conditions. 

After  establishing  a  single-variable  non-linear  method  to  produce  the  initial 
conditions,  the  process  was  run  with  various  chief  inclinations.  The  results  showed  that 
the  method  was  consistent  at  producing  secular  drifts  of  approximately  two  to  five 
centimeters  drift  in  the  orbital  direction  over  the  25  orbit  propagation  period. 
Unfortunately,  the  secular  drift  rate  in  the  out-of-plane  direction  was  less  consistent 
across  the  inclinations,  with  the  worst  drift  being  1.4  centimeters  per  orbit  at  an 
inclination  of  30°  and  the  best  being  approximate  41  micrometers  per  orbit  at  95°.  This 
suggests  that  the  ability  of  this  method  to  reduce  the  secular  drift  in  the  z-direction  is 
somewhat  dependent  on  the  inclination.  According  to  the  paper  done  by  Sabatini  (5:97- 
100)  this  is  also  the  case  with  the  h  invariant  method  developed  by  Schaub  and  Alfriend 
(6:77-95).  However,  the  previous  method  showed  most  desirable  results  at  the  critical 
inclination,  whereas  this  method  produces  optimum  results  at  near  polar  orbits. 
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When  the  method  was  extended  to  include  non-circular  chief  orbits,  the  method 
was  less  successful  at  producing  small  drifts  in  the  out-of-plane  direction;  producing  drift 
rates  of  over  3  centimeters  per  orbit.  However,  it  was  shown  that  the  method  was  still 
able  to  produce  relative  stability  in  the  x-direction,  with  less  than  2  centimeters  drift  after 
propagating  the  initial  conditions  for  25  orbits.  Unfortunately,  when  attempting  to 
control  the  initial  relative  angle  [>  the  size  of  the  relative  orbit  increased  by  over  100  km. 
This  is  due  to  the  moderately  large  eccentricity  of  the  chiefs  orbit  and  the  difference  in 
the  argument  of  perigee  required  to  induce  the  initial  relative  angle.  This  suggests  that  in 
order  to  maintain  flexibility  of  the  deputy’s  position  on  relative  orbit  the  chief  will  have 
to  be  at  a  near  zero  eccentricity. 

After  determining  that  the  single-variable  method  was  unable  to  produce  desirable 
drift  rates  for  all  inclinations,  the  single -variable  method  was  then  extended  to  solve  for 
both  variables  simultaneously.  This  two-variable  non-linear  method  produced 
significantly  more  accurate  results  for  the  lower  inclinations;  however,  even  with 
increased  accuracy  in  the  lower  inclinations,  the  method  was  not  always  able  to 
reproduce  the  accuracy  established  by  the  single-variable  method  for  near  polar  orbits. 
This  suggests  that  if  it  is  possible  to  place  a  formation  at  inclinations  around  or  above  the 
critical  inclination,  one  of  these  two  methods  will  be  able  to  establish  near  h  invariance. 
Fortunately,  even  if  this  is  not  possible,  these  results  do  show  that  the  two-variable  non¬ 
linear  method  will  be  able  to  significantly  reduce  the  secular  drifts  due  to  the  perturbing 
accelerations. 

Future  Work 

Further  work  on  this  subject  could  include  considering  the  stability  of  multiple 

satellite  formations.  By  declaring  certain  deputies  as  “sub-chiefs”  it  would  be  possible  to 
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determine  the  ROE  between  the  sub-chiefs  and  the  deputies.  This  could  be  repeated  until 
the  relative  geometry  between  all  of  the  satellites  in  the  formation  has  been  detennined. 
This  would  allow  for  a  better  understanding  of  more  advanced  formations;  where  the 
relative  geometries  between  the  chief  and  deputies  are  a  concern,  as  well  as  the  relative 
geometries  between  deputies.  This  information  could  also  simplify  the  task  of 
deconflicting  the  trajectories  of  multiple  deputies,  and  decreasing  the  possibility  of  a 
collision.  Additional  work  could  attempt  to  further  refine  the  relative  positions  in  an 
attempt  to  achieve  specific  geometric  shapes  at  a  desired  argument  of  latitude.  This 
could  aid  in  spreading  the  system  out  normal  to  a  point  of  interest,  and  allow  for  a 
broader  distributed  aperture  system. 
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Appendix  A:  State  Transition  Matrix 
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Appendix  B:  Figure  of  Relative  Orbital  Elements 


Definition  of  Relative  Orbital  Elements 
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Figure  29:  Figure  of  Relative  Orbital  Elements 
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